Graphene on Ir(lll): Physisorption with chemical modulation 
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The non-local van der Waals density functional (vdW-DF) approach is applied to calculate the 
binding of graphene to Ir(lll). The precise agreement of the calculated mean height h = 3.41 A 
of the C atoms with their mean height h = (3.38 =b 0.04) A as measured by the X-ray standing 
wave (XSW) technique provides a benchmark for the applicability of the non-local functional. We 
find bonding of graphene to Ir(lll) to be due to the van der Waals interaction with an antibonding 
average contribution from chemical interaction. Despite its globally repulsive character, in certain 
areas of the large graphene moire unit cell charge accumulation between Ir substrate and graphene 
C atoms is observed, signaling a weak covalent bond formation. 

PACS numbers: 68.43.Bc, 68.65.Pq, 68.49.Uv, 73.20.Hb 



Epitaxial growth on metals is a key method to pro- 
duce high quality graphene on large scales [TJ [2] . Owing 
to the strength of the C-C bonds, large (incommensu- 
rate or weakly commensurate) superstructures are found 
for lattice-mismatched systems. The extend of superpe- 
riodicity in the electronic structure [3j [4] and buckling of 
the carbon layer [5j [6] depends on the strength and local 
variation of the interaction between graphene and sub- 
strate, which ranges from strong chemisorption to weak 
physisorption |7J. A key parameter for this strength is 
the height h of the carbon adsorbate, which is, however, 
difficult to assess by experiment: Analysis by x-ray or 
electron diffraction necessitates a large number of fitting 
parameters because of extended unit cells [8 , 9J. In addi- 
tion, the x-ray scattering amplitude of C is low [9J. Deter- 
mination of h by scanning tunneling microscopy (STM) 
and atomic force microscopy is also questionable [10J. 
Only for a few systems such measurements have been 
performed, and consensus has not been reached [see [6] 
for Ru(0001)], except for highly commensurate graphene. 

In density functional theory (DFT) a quantitative 
description of the interaction between graphene and 
metal is a challenge, because the most commonly used 
exchange-correlation functionals [local density approxi- 
mation (LDA) and generalized gradient approximation 
(GGA)] are (semi)local [11]. They lack the nonlocal- 
correlation effects responsible for van der Waals (vdW) 
interaction [12 , 13J. As an example, for the weakly bound 
system graphene/Ir(lll) investigated here, DFT-GGA 
calculations result in a mean value of h ~ 3.9 A [5 J 
and very low binding energies of only a few meV 
per C atom [5j [14], thus contradicting the experimen- 
tally observed formation and stability of graphene at 
about 1300 K [lj. Such a low E^ is also not plausible 



in view of the observation that the peeling force needed 
to remove graphene from Ir(lll) [15J is higher than what 
is necessary to exfoliate graphite (experimentally deter- 
mined E\> — — 52meV per C atom [H]). LDA calcula- 
tions [T71 [T8] succeed in binding graphene in accord with 
h = 3.42 A. However, conceptually this success is un- 
satisfying and questionable as the error of disregarding 
the vdW interaction is diminished by error cancellation: 
The LDA has a well known tendency to overbind. Only 
recently, the first truly nonlocal correlation functional 
vdW-DF [3] was developed and successfully applied to 
simple vdW-bonded systems [20J , thereby opening a new 
perspective to correctly describe also complex systems 
with significant vdW bonding. 

Without a proper understanding of binding between 
graphene and Ir(lll), also other discrepancies will be 
hard to resolve: Whereas previously the absence of any 
interaction between the Dirac cone and the Ir 5d bands 
was assumed [3j, more recent experiments point towards 
a hybridization close to the Fermi level [4J. Such a hy- 
bridization also influences the properties of phonons in 
graphene: Strikingly, both the presence [21] and the ab- 
sence [4J of the respective Raman modes have been re- 
ported. Even the commonly assumed shape of graphene 
[5] appears inverted in a recent AFM experiment [22J. 

In this Letter we further the understanding of graphene 
on metal by (i) performing an x-ray standing wave 
(XSW) experiment, which allows us to unambiguously 
specify the average bond distance h and the amplitude 
of buckling Ah and (ii) applying nonlocal vdW-DF using 
a realistically large supercell. The comparison of exper- 
iment and theory allows to benchmark the functionals 
used to describe bonding of graphene and the bonding of 
7r-conjugated systems on metals in general. 
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The experiments were performed at the ESRF un- 
dulator beamline ID32 [23]. Ir(lll) (mosaicity 0.2° 
determined from x-ray diffraction rocking curves) was 
cleaned in UHV by cycles of 1.5 keV Ar + bombardment 
and annealing to 1420 K. Graphene was grown by re- 
peated cycles of room temperature C2H4 adsorption fol- 
lowed by thermal decomposition at 1420 K (temperature 
programmed growth or TPG). Well oriented graphene 
forms, leading to the characteristic moire pattern in low 
energy electron diffraction [3]. We prepared a cover- 
age of (0.39 ± 0.03) ML via two TPG cycles and of 
(0.63 ±0.04) ML via four cycles. For coverage estimation 
we used the fact that a fraction of (0.22 ± 0.02) of the 
free Ir(lll) surface is covered with graphene after each 
TPG cycle [24]. Photoelectron spectra (overall resolution 
400 meV) were recorded using a hemispherical electron 
analyzer. The reflectivity of the sample was measured by 
directing the reflected beam to an insulated metal plate 
and measuring the resulting electron emission current. 
Scanning tunneling microscopy (STM) was performed on 
another Ir(lll) sample in a separate UHV system where 
the growth was repeated [25J. 

An X-ray standing wave was created in the interface 
region of a crystal using Bragg reflection. We used the 
(111) reflection at an angle close to 90° at 2.801 keV. The 
XSW maxima, which are periodic with the Ir(lll) lat- 
tice planes, are shifted by half the lattice plane distance 
when scanning through the Bragg reflection by changing 
the beam energy. The photoelectron yield of adsorbates 
during such a scan depends on their height above the sur- 
face. From the results of an XSW measurement, one ob- 
tains the coherent position P H and the coherent fraction 
F H [3j, which roughly correspond to the mean adsorbate 
height and the spread around this value [25J. 

From the photoelectron spectra [compare Fig.[TJa)] we 
determined the Cls binding energy to (284.2 ± 0.1) eV, 
in agreement with Ref. 7. After subtraction of a Shirley 
background, all Cls spectra could be fitted well with a 
single Gaussian. From the resulting peak area as a func- 
tion of photon energy [Fig. [ljb)] and averaging over sev- 
eral scans for the same preparation, the structural pa- 
rameters P H = 0.53 ± 0.01 and F H = 0.87 ± 0.04 for 
0.39 ML and P H = 0.52 db 0.01 and F H = 0.74 ± 0.04 
for 0.63 ML are obtained [23J. To interpret these val- 
ues, we tested three simple height distribution functions 
(Gaussian, rectangular, p6m layer with Fourier compo- 
nents up to first order [25J). It turns out that almost 
the same structural parameters give the best fits to P H 
and F H for all three models, showing the robustness of 
our interpretation. We determined ft = (3.38 ± 0.04) A 
[^5] for both coverages, and a standard deviation of 
<r h = (0.19±0.03) A for 0.39 ML and a h = (0.27±0.04) A 
for 0.63 ML. Here the errors contain both the experi- 
mental uncertainty as well as the small deviations result- 
ing from the choice of the model. The standard devi- 
ation gives an upper limit for the possible corrugation 
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FIG. 1: (Color online) (a) Squares: X-ray photoelectron in- 
tensity of the Cls-peak for 0.63 ML of graphene on Ir(lll). 
Thin line: Shirley-type background, Thick line: Gaussian fit 
(full width at half maximum 0.85 eV) above background, (b) 
X-ray reflectivity (triangles) normalized using the fit to the 
data (solid line [27] ) and exemplary Cls photoemission yield 
normalized to an off-Bragg yield of 1 (squares: 0.39 ML, cir- 
cles: 0.63 ML, shifted upwards by unit of 1 for easier readabil- 
ity). Solid lines: fits to the data [231 I27| taking into account 
additional broadening due to the mosaicity of the sample and 
the bandpass of the X-ray beam. 

of graphene. The measured mean height is similar to 
the interlayer distance in graphite of 3.36 A. This al- 
ready indicates that only weak bonding can be at place. 
Where applicable, the standard deviation translates to 
Aft = (0.6 ± 0.1) A for 0.39 ML and Aft = (1.0 ± 0.2) A 
for 0.63 ML, see discussion below. 

Ab initio DFT calculations have been performed by 
using the projector augmented wave method [28J with 
the PBE-GGA functional [29J as implemented in VASP 
[5U] . To obtain a reliably relaxed adsorbate-surface geom- 
etry [31] we have implemented vdW forces using the semi- 
empirical method DFT-D [32J , where the Ir Cq coefficient 
was determined by comparing DFT-D and nonlocal vdW- 
DF calculations for several adsorption geometries of ben- 
zene on Ir(lll). For the final relaxed graphene-Ir(lll) 
geometry, the total energy was calculated with vdW-DF 
in a post-processing approach using the JuNoLo code 
[33J. In this method, the vdW-DF total energy is eval- 
uated with the electron density obtained from the GGA 
calculations and its value is not changed when performing 
a full selfconsistent vdW-DF calculation [34] . 

The system was modeled by (10 x 10) unit cells of 
graphene on (9 x 9) cells of Ir(lll) resulting in a ratio 
of the surface unit cell lengths of ac/a\ T = 1.111, close 
to the experimental value of acja\ x = 1.107 for flakes 
with an average size of 1000 A. We used a slab geometry 
including in total about 600 atoms (4 layers of Ir and 1 
layer graphene) in the unit cell. In previous studies for 
graphene on metals the vdW-DF functional was applied 
using much smaller unit cells [T2| \13\ [35] where effects 
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FIG. 2: (Color online) (a) Top view and (b) side view [cut along the dashed line in (a)] of the relaxed structure of 
graphene/Ir(lll) obtained by DFT including vdW interactions. Regions of high-symmetry stacking (fee, hep, top) are marked 
by circles (a) or arrows (b-d). The color scale in (b) and (a) ranges from h = 3.20 A (dark) to h = 3.65 A (yellow), (c) 
Visualization of the nonlocal correlation binding-energy density ec l,hind (r) caused by adsorption. The color scale ranges from 
e^' bind (r) = meV A~ 3 (blue) to e? 1,bind (r) = -28.3 meV A~ 3 (red), (d) Charge transfer upon adsorption. The color scale 
ranges from Ap = —0.0138 e A -3 (blue) to Ap = 0.0138 e A -3 (red). A negative value indicates loss of electron density, (e) 
Magnified view of red box in (d) Views from different angles can be obtained in [25] . 



of the moire superstructure are suppressed. Plane waves 
with a kinetic energy up to 400 eV have been included 
in the basis set and the Brillouin zone was sampled by a 
3x3x1 /c-mesh. Geometry optimization was obtained 
by relaxing the top 2 Ir layers and the graphene layer 
including the long-range vdW forces as described above 
with a force threshold set to 1 meV/A. 

The geometry resulting from this calculation is shown 
in Fig. ^[a). The largest height of 3.62 A is found in re- 
gions where the center of the C hexagon is located on top 
of an Ir atom (top region) and the lowest height of 3.27 A 
in the hep-region (center of C hexagon above threefold co- 
ordinated hep site), slightly lower than the one of 3.29 A 
in the fec-region. The mean height is h = 3.41 A, in excel- 
lent agreement with the experimental value. The corru- 
gation resulting from DFT including vdW is Ah = 0.35 A 
or = 0.09 A, and thus safely within the upper limits 
determined by experiment for both coverages. 

The averaged binding energy per C atom in our cal- 
culations is E h = —50 meV/C. The vdW-DF approach 
makes it possible to distinguish between local and non- 
local contributions to the overall binding energy [25J . The 
nonlocal correlation energy E£ l can be expressed as a 
function of the charge density p(r) 

Ef = \JJ drdr'p(r)<j>(r,r')p(r') = J dref{r), 

where 0(r, r') is the kernel function (see discussion in 
Ref.S|) and e^(r) is the nonlocal correlation energy den- 
sity at each point r in real space. 

The distribution of the nonlocal correlation binding en- 
ergy density e£ 1,bmd (r) (i.e. the change in e^(r) caused 
by adsorption [25]) is shown in Fig. [2] (c). Note that 
e" 1 ' bmd (r) arises from a non-local quantity as its value 



at a given point depends on the interaction of the charge 
density at this point with the one at all other points. 
The fact that this energy density is broadly distributed 
in a layer just above the metal surface and just below the 
graphene plane (note that for graphene we cut alternately 
through bonds and hexagon centers) shows that the po- 
larization effects responsible for the vdW interaction are 
spread over the entire sheet, which is a clear fingerprint 
of a vdW-bonded 7r-conjugated system [T^l I5T]. 

For the relaxed geometry shown in Fig. [5|a), the 
binding energy E^ GA calculated in GGA is repulsive 
+20 meV/C) while the non-local binding energy con- 
tribution E£ l « — 70 meV/C is attractive summing up to 
the total vdW-DF binding energy E h = -50 meV/C [25]. 
Contrary to the impression given by these averaged val- 
ues, the binding is not pure physisorption, but chemi- 
cally modulated. This becomes obvious when analyzing 
the charge transfer caused by adsorption [Figs. |2|d) and 
(e)]. In the hep- and fec-regions a small charge transfer 
from graphene towards the substrate takes place. A C 
atom sitting directly atop of an Ir atom [see Fig. [2|e)] 
hybridizes its C(2p 2 ) orbital with the Ir(5d 3 ^2_ r 2) or- 
bital. As a result charge accumulates just in between the 
C atom and the Ir atom, indicating formation of a weak 
covalent bond. This charge is provided primarily from 
the neighbor C atoms to the bondforming ones sitting 
atop of Ir atoms. The charge deficit of these neighbors 
explains their tendency to bind additionally deposited 
metal atoms [5J. This charge transfer is intimately re- 
lated to the nonlocal binding energy density [Fig. [2|c)] 
localized in specific regions close to those Ir atoms where 
the charge transfer from graphene to substrate occurs. 
This indicates those sites at the metal substrate that 
become more polarizable upon adsorption due to the 
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graphene-surface interaction. In total, graphene has lost 
« 0.01 electrons/C of charge resulting in slight p-doping. 

A by-product of DFT calculations are the Kohn-Sham 
eigenstates. Due to the large supercell in our calculations 
the calculated bands are multiply folded. To extract the 
dispersion relations for the entire Brillouin zone heavy 
postprocessing would be necessary [36]. Nevertheless, 
taking only the states with a large projection on the car- 
bon atoms, thereby filtering out the substrate states, we 
identify the Dirac cone of the adsorbed graphene in the 
vicinity of the Fermi level. The Dirac point is shifted 
0.2 eV above 

-^Fermi consistent with the experimental 
value of 0.1 eV [3J. 




FIG. 3: (Color online) STM topographs of graphene on 
Ir(lll). (a) 0.39 ML. Flakes with the characteristic graphene 
moire are visible, (b) 0.63 ML. A sheet extends over a wide 
area including a substrate step. Lowest regions are bare 
Ir(lll). Image widths 1000 A, bias voltage [/ samp ie = — 1 V, 
tunneling current I = 0.2 nA (a), I = 0.06 nA (b). 

Finally we address the experimentally observed cov- 
erage dependence of the corrugation. To obtain insight 
into this dependence we analyzed the graphene films us- 
ing STM (Fig.[||). For 0.39 ML (Aft = 0.6 A), graphene is 
present in the form of flakes with a mean size of ~ 500 A, 
most of them located on extended terraces. For 0.63 ML 
(Aft = 1.0 A) large coalesced flake agglomerates form 
with linear dimensions above 1000 A. Here, the graphene 
overgrows substrate steps [lj and occasionally wrinkles 
are found [37J, which are absent for the smaller flakes. We 
speculate that in the different geometries the shrinking 
of the substrate while cooling down from the high growth 
temperatures [37J has distinct effects: As small flakes are 
able to float they remain relaxed and flat. Larger flakes 
pinned to and overgrowing steps are unable to float. One 
mechanism to avoid build up of stress is then to buckle. 
Dedicated experiments to substantiate these speculations 
are necessary and under way. Such a strain dependence 
of the corrugation is out of reach in our simulations, as 
in the vdW-DF calculations the 10 x 10 graphene sheet 
matches the 9x9 substrate mesh without strain. 
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GRAPHENE ON IR(lll): PHYSISORPTION VERSUS CHEMICAL REPULSION - SUPPLEMENTARY 

MATERIALS 

Scanning tunneling microscopy 

The STM measurements depicted in Fig. 3 of our manuscript have been performed in a UHV system at the 
University of Cologne using a home-build beetle-type STM (the system is described in Detail in Refs. [TJ |2j). For the 
measurements we used an Ir(lll) crystal from the same manufacturer with identical specifications and prepared the 
crystal and the graphene using the same experimental parameters as for the XSW experiments. 



Analysis of the X-ray standing waves data 



Since we detected the differential photoelectron yield, multipole contributions had to be taken into account and the 
data was fitted according to |3J . Strictly, XSW can only give you the relative height with respect to the crystal lattice 
planes, i.e. h mod di r (iii)- Here we disregarded unphysically small or large values. To interpret the experimentally 
determined parameters P H and F H , we tested different height distribution functions G(z): A Gaussian function: 



G(z) 



1 



(i) 



a rectangular function with width s and height 1/s (note that for such a distribution the standard deviation is 
given as a = J j^s) 



G(z) 



(2) 



- s if|*-*b|<§, 
if|z-z |>§, 

and the height distribution of a sixfold symmetric (p6m) layer described by a superposition of three cosine waves 
(eggbox- model), i.e. containing only Fourier coefficients up to first order: 



2 / 

h(f) = h + ^cos(fcir) + cos(A?2r) + cos(A?3r) 



(3) 



with |fci| = \k 2 \ = |fc 3 | and lk u k 2 = lk 2 ,k 3 = 120°, 

The eggbox-model is motivated by the appearance of C/Ir(lll) in STM. Note that the exact shape resulting from 
the calculations can well be fitted with this eggbox-model with a root mean square deviation of only 0.006 A/C. When 
in the studies for graphene on Ru(0001) cited in the main manuscript only the highest and lowest value of h were 
explicitly given, we estimated h on the basis of the eggbox-model. 



Details of the calculation 



The vdW-DF approach makes it possible to distinguish between local and non-local contributions to the overall 
binding energy. Technically, the correlation part of the generalized gradient approximation (GGA) functional is 
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replaced by a sum of the local density approximation (LDA) correlation functional £^lda,c an d a nonlocal term E^ 1 
such that the total energy is given by 



^vdW-DF = ^GGA — ^GGA,c + ^LDA,c + E£ l , 



(4) 



where Eqqa is the self-consistent GGA total energy evaluated in a conventional DFT calculation. The nonlocal 
correlation energy E£ l can be expressed as a function of the charge density n(r), 



Ef = \J J drdr'n(r)^(r,r')n(r'), 



where the kernel function 0(r, r') is discussed in detail in Ref. [4j 
To evaluate Eq. [5| it can be rewritten as 



(5) 



E. 



nl 



drdr'n(r)<p(r, r')n(r'), 
= J drn(r) ^ J dr'<j>(r, r')n(r') 
= J drn(r)ef{r), 
dref(r), 



where e^ l (r) is the nonlocal correlation energy density for each point r in the real space. 

Then, the nonlocal correlation binding energy density e£ 1 ' blnd (r) plotted in Fig. 3 is defined as 



(6) 



e^' bind (r) = e^' sys (r) - ( e ^ a P heM ( r ) + e ni,ir(m)( r )) 5 ( 7 ) 

where the e^' sys (r) denotes the nonlocal correlation energy density of the graphene-Ir(lll) system, while 
e ni,graphene^ r j an( j e £ 1,Ir ( m ) ( r ) are the nonlocal correlation energy density at point r of the graphene and Ir(lll) in 
their graphene-surface relaxed geometry, respectively. The non-local binding energy contribution, E^ 1 , per C atom, 
is given by the integal of e£ 1,blnd over the space of the unit cell, plus a small contribution arising from the binding 
energy differences of ^gga,c + ^lda,c in Eq. [4j divided by the number of C in the unit cell. 
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